Updates in SASfit for fitting analytical expressions and numerical models to small-angle scattering patterns

Recent enhancements and additions to the SASfit program are discussed, including anisotropic scattering models, flexible distributions, regularization techniques such as the expectation-maximization method, and new structure factors, especially for ordered nano- and meso-scaled material. The Ornstein–Zernike solver for numerical structure factors is also introduced.


OPO -oriented primitive objects
For many different geometric objects such as sphere, toroid, cone, pyramid or tetrahedron, analytical form factors are made available in SASfit. For analyzing anisotropic scattering patterns from such geometries SASfit expects one-dimensional SAS profiles in radial direction at a certain angle or in azimuthal direction.
A subclass of form factors in SASfit are oriented primitive objects. At the beginning this group was intended for simple single phase objects like homogeneous cube, cylinder, ellipsoid and their affine deformed shapes. We than added some other objects like platonic solids and super quadrics and further extended it by cones, pyramids, etc. All of them can be either perfectly oriented in the beam or are fully randomly oriented.
For the perfect oriented version three Euler angles need to be supplied. Also many input parameters are needed for describing the affine deformation. The most simple geometric objects are unit spheres, unit cube and unit cylinder, and their affine deformed shapes of ellipsoid, oblique parallelepiped and oblique cylinder for which their Q-dependent form factors F (Q) are given analytically.

Ellipsoid
Parallelepiped Cylinder where η(r) is the scattering length density distribution of the scattering object and V the sample volume. In the following it is assumed that the particles have a homogeneous scattering length density η P , embedded in a matrix of also homogeneous scattering length density η M = η P − ∆η. The position, orientation and size of the objects is uniquely defined by the center of scattering length density R, the half axes det(D −1 ) = abc (e ax e by e cz + e ay e bz e cx + e az e bx e cy − e az e by e cx − e ay e bx e cz − e ax e bz e cy ) (7) whereas det(D) det(D −1 ) = det(D D −1 ) = det(1 I) = 1 To reorient or rotate the objects the direction of the base vectors e a , e b , e c needs to be adjusted which is done easily by introducing a rotation matrix R αβγ . The rotation matrix can be defined via Euler angles α, β and γ. Its definition depends on the convention used. There exist twelve possible sequences of rotation axes, divided in two groups: proper Euler angles (z-x-z, x-y-x, y-z-y, z-y-z, x-z-x, y-x-y) and Tait-Bryan angles (x-y-z, y-z-x, z-x-y, x-z-y, z-y-x, y-x-z). The difference between the two conventions is that Tait-Bryan angles represent rotations about three distinct axes, while proper Euler angles use the same axis for both the first and third elemental rotations. SASfit allows internally to use any of the 12 conventions. However, as default convention the yaw-pitch-roll is used (also named gier-nick-roll, east-north-up, north-East-Down or IUCr macros version 2.1.15: 2021/03/05

Randomly oriented objects
Next to the form factor of an oriented primitive object also the form factor of a fully randomly oriented version is made available. The orientational average is done by integrating over the full Ewald's sphere, which is a double integration. For this orientational average SASfit provides optimized integration routines which can be configured in the GUI shown in figure 1 of the article. Next to the multidimensional adaptive integration routines pcubature and hcubature (Johnson, 2020) also some other routines are supplied. However, these routines are nonadaptive without an error estimate of the integral but optimised for integrations over a sphere surface. Nevertheless, these routines often need significant less function evaluations than conventional multidimensional integration routines. The implemented methods are 2D-GaussLegendre integration, Lebedev (Lebedev, 1975;Lebedev, 1976;Lebedev, 1977) and Fibonacci (Niederreiter, 1992;Marques et al., 2013) quadrature algorithm with a fixed and configurable number of function evaluations and without error estimation of the integral.

Size distributions
A standard task in analysing small-angle scattering data is the determination of the pair-distance distribution function or a size distribution of particles with a known form factor. Also desmearing of small-angle scattering data belongs to the same family of problems. The size distribution can be determined if the particle shape is known. For many particle shapes the scattering intensity or form factor can be calculated via an analytical expression. The scattering intensity is then described by integrating the squared form factor F 2 (Q, s) over a distribution function N (s), where s describes the size of the scatterers (Pedersen, 1997): To solve equation 25 often the size distribution is assumed to be an analytical distribution function such as LogNormal, Gaussian, Weibull or Schulz-Zimm with few input parameters. These parameters are then optimised by performing the numerical integration and applying a Levenberg-Marquardt strategy to solve a nonlinear least squares problem (Breßler et al., 2015). This however assumes a priori knowledge about the distribution function.

MetaLog distribution
A distribution function designed to be very flexible is the metalog distribution. The quantile function of the metalog distribution Q(y) = M k (y) is defined as (Keelin & Powley, 2011;Keelin, 2016;Keelin et al., 2019) The metalog PDF as defined above is an unbounded distribution. In case of describing a size distribution it needs to be bounded (Logit metalog) or at least semibounded (bounded below, Log metalog). This can be easily done by a transformation a lower bound b l and an upper bound b u , where z is then metalog-distributed. The corresponding quantile functions and PDF read as

Regularization
Instead of using an analytical description of the size distribution it would be more desirable to use a discretization of the size distribution N i at the positions s i and determine all N i without any further assumptions.
Mathematically, these problems are equivalent to finding a solution for a Fredholm integral of the first kind. These tasks are ill-posed problems by their nature and an additional regularization technique has to be included in the analysis to obtain a stable solution (Tikhonov, 1943;Tikhonov et al., 1995). For small-angle scattering this technique has been introduced by (Glatter, 1977) and later on taken over by (Svergun et al., 1988) and (Hansen & Pedersen, 1991). To find a stable solution a penalty term in the minimization procedure had to be included which introduced the question how to properly weight this penalty term. An early discussion about this for small-angle scattering can be found in (Svergun, 1992). However, the proper weighting of the penalty has been treated also in other fields and a discussion can be found for example in (Hansen & O'Leary, 1993;Hansen, 2000). Fredholm integrals, like the one in eq. 25, can be written as linear matrix equations Ax = b by discretization. Hereby, A is a M × N matrix, b is a vector of length M and x is a solution vector of length N . Their elements are given by The solution vector x is then obtained by a least squares optimization. If the penalty term L can also be linearized the problem can be reduced to solving the linear least squares problem where || . . . || 2 w denotes the weighted sum of squared residuals and λ the regularization parameter weighting the penalty term. Typically the linear least squares equation has to be solved for many λ and a criteria of finding the proper regularization parameter has to be chosen (Hansen & O'Leary, 1993). For the linear penalty term typically the identity matrix or the first or second order derivative operator are taken (Donatelli & Reichel, 2014). Using the entropy with a positive prior as a penalty function the solution vector is also automatically constrained to positive values. Maximum Entropy methods have been used for small-angle scattering analysis, for example by (Hansen & Müller, 1996;Hansen & Pedersen, 1991). (Steenstrup & Hansen, 1994) have shown how to use Maximum Entropy methods without the positivity constraint.

Expectation Maximization
Expectation maximization (EM) algorithm is a procedure for designing iterative methods for likelihood maximization. This framework has been explained first by (Dempster et al., 1977). Applied to Fredholm integrals the method is an iterative fixed point method for non-negative real valued functions (Vardi & Lee, 1993). The method described there is equivalent to the Lucy-Richardson method (Richardson, 1972;Lucy, 1974) which is well known and widely used in image analysis.
For inverting scattering data to pair distance distribution functions the condition of non-negativity of the involved functions does not hold anymore as both the kernel of the Fredholm integral as well as the pair distance distribution function can take negative values. In (Chae et al., 2018) it has been shown, how the EM algorithm can be reformulated for non-density functions, i.e. to functions which also can take negative values. The EM iteration scheme or Lucy-Richardson inversion method reads as where O EM is the fixed point operator. It has been shown, that the convergence is assured since the algorithm is guaranteed to increase the likelihood at each iteration.
On the other side the algorithm is known to converge quite slow.  test, which is obtained by a second order Taylor expansion of the natural logarithm around 1 of the G-test.
weighted sum of squared residuals (38) with I exp (q i ) = I exp (q i )/ i I exp (q i ) and I th (q i ) = I th (q i )/ i I th (q i ). Furthermore, the roughness of the resulting size distribution is plotted which has been described by either the sum of the squared first or second derivatives of all radii, i.e. i |dN  Smoothing Operator In SASfit the smoothing operation suggested by (Eggermont, 1999;Eggermont & LaRiccia, 1995;Byrne & Eggermont, 2011) has been implemented.
They suggest to add before or/and after each normal EM step (eq. 36) an additional smoothing step so that the EM algorithm takes the form.
with the smoothing operator whereby 0 < h 0.3. In case of a double smoothing they suggest a nonlinear smoothing operator before each EM step so that a complete smoothed EM step reads as    Figure 7 shows the logarithm of a roughness term, namely the sum of 1 st and 2 nd derivatives of the resulting size distribution and the entropy against the goodness-of-fit parameter from eq. 40, namely the "G-test" parameter. The optimum smoothing parameter has been defined via the corner of the L-curve defined by (log(G), log(||Lx|| 2 )), withL being the first or second-order discrete derivative operator. for example given in (Donatelli & Reichel, 2014) and by (log(G), log(S)) with S being the entropy defined in eq. 47. The corners were determined numerically as they are defined by having the maximum curvature. They are marked in fig. 10 by a symbol. The whole analysis does not make use of the uncertainties of the measured intensities. If the error bars are trustable and the model for the form factor is correct the weighted sum of squared residuals should be χ 2 r = 1.
This can also be used as a condition for determining the smoothing parameter h. In the case of the present analysed data both strategies yield similar results as can be seen in fig. 11.
Penalty function The Lucy-Richardson method (Richardson, 1972;Lucy, 1974) has been extended for an additional penalty function, namely the entropy (Lucy, 1994). The maximum entropy penalty can be introduced in the iteration algorithm by knowing either a fixed or assuming an adaptive prior for the solution vector. The entropy S is given by with m being the prior estimate of x. The update scheme in eq. 36 has to be extended as in case of a known prior. If a prior is not known a constant prior can be chosen. However, even if no model for the prior is known, one can try to construct one adaptively according to Horne (Horne, 1985) by applying for example a Gaussian point spread function to the actual solution and using this as a prior for the next iteration step. In this case the update scheme changes to with a smoothing operator Π ij for the prior:  fig. 13 together with the exact solution and for a Lagrange parameter leading to χ 2 r = 1. The L-curves in fig. 12 show that the G-test parameter for the constant prior varies over several orders of magnitudes whereas the adaptive prior, even with a Lagrange parameter far away from its optimal value leads to a reasonable good result. Therefore, the L-curve shows a much smaller range on the x-axis for the adaptive prior.  fig. 12 and the solution of the EM iteration scheme using Maximum Entropy term with constant prior (left) and adaptive prior (right) is compared with the exact solution as well as for a Lagrange parameter leading to a weighted sum of squared residuals with χ 2 r = 1.
Finding fixed points efficiently There are several very efficient algorithms for increasing the convergence rate of contractive fixed point operators: Lewitts method Lewitt et al. (Lewitt & Muehllehner, 1986) have suggested to increase the convergence rate by introducing an over-relaxation parameter in the iter-ation loop: where λ (k) is the over-relaxation parameter (λ (k) > 1) whose purpose is to accelerate the iterative process, and whose value is not so large as to violate the non-negativity condition x (k+1) j > 0. To ensure that the over-relaxation of negative corrections to the solution vector does not decrease the value of any solution vector component below zero they suggest the following algorithm: first a pseudo-relaxation parameter µ Then they defined a critical relaxation parameter as follows: Finally the over-relaxation parameter λ (k) was chosen according to with λ max 4 so that 1 ≤ λ (k) ≤ λ max .
Anderson acceleration Another more efficient method to accelerate the convergence rate of the EM iteration scheme is the so called Anderson acceleration (Anderson, 1965) for solving fixed point problems [see also (Walker & Ni, 2011;Toth & Kelley, 2015)], it has also been suggested for the EM iteration scheme (Henderson & Varadhan, 2019). In SASfit several fixed point accelerations are supplied by making use of the sundials library (Hindmarsh et al., 2005). However, only the Anderson acceleration speeds up the fixed point iterations in a reliable way whereas all the others such as GMRES (Saad & Schultz, 1986), FGMRES (Saad, 1993), Bi-CGStab (van der Vorst, 1992) or TFQMR (Freund, 1993) become unstable and do not converge.

Structure factors from solving Ornstein-Zernike equations
The central part of the SASfit algorithm to obtain a structure factor from a given pair interaction potential performs an iterative solution of the Ornstein-Zernike equation (Ornstein & Zernike, 1914). The theoretical background can be found in great detail in (Nägele, 2004;Borowko et al., 2000;Hansen & McDonald, 2013;Bomont, 2008;Caccamo, 1996;Likos, 2001). The Ornstein-Zernike equation reads h(r 12 ) = c(r 12 ) + c(r 13 )ρ(r 3 )h(r 32 )dr 3 + · · · where h(r 12 ) is the total correlation function, c(r 12 ) the direct correlation function (direct effect of particle 1 on particle 2), c(r 13 ) the direct correlation function describing the effect of particle 1 on particle 3 which influences particle 2, and ρ the particle number density of the colloids, molecules, or atoms that form the liquid. If the fluid is uniform and isotropic, the Ornstein-Zernike relation becomes The basic idea of this equation is that the total correlation between positions of two particles is a combination of their direct and indirect (through neighboring particles) interactions. The first contribution to h(r) is the direct correlation function c(r) that represents the correlation between a particle of a pair with its closest neighbor separated by a distance r. The second contribution is the indirect correlation function γ(r), which represents the correlation between the selected particle of the pair with the rest of the fluid constituents. g(r) is known as the radial or pair distribution function, which measures the probability that given a particle at the origin, another particle of the fluid can be found at a distance r from it. When the distance separating a pair of particles tends to infinity, the correlations vanish and g(r) tends to 1. This means that the total correlation function defined as h(r) = g(r) − 1 tends to 0. The Fourier transform F{} of g(r) and h(r) are directly related to the structure factor S(q) that is experimentally measurable by x-ray or neutron scattering and S(q) = 1 + ρh(q) whereh(q) = F {h(r)} is the Fourier transform of the total correlation function h(r), so thath Via inverse Fourier transform F −1 {} of the structure factor the pair correlation function is obtained As for g(r) also the direct correlation function c(r) can be, in principle, derived from the experiment as By Fourier transformation of eq. 61, one obtains wherec(q) = F {c(r)} is the Fourier Transform of c(r). In order to determine the two correlation functions h(r) and c(r) for a given pair potential u(r), eq. 61 must be supplemented by an auxiliary closure relation. Furthermore the Fourier transform of eq. 62 in combination with eq. 69 allows to write the Fourier transform of the indirect correlation function as The closure is not complete as long the bridge function is unknown. Even though the bridge function is exactly defined the function is an unknown function of inter-particle distance and only approximations can be given. The direct correlation function can be written in terms of the indirect correlation function, the interaction potential and the bridge function as and h(r) + 1 = g(r) = exp [−βu(r) + γ(r) + B(r)] where u(r) is the pair interaction potential, β = 1/(k B T ) the inverse temperature, with k B being Boltzmann's constant, T is the absolute temperature, γ(r) being the indirect correlation function and ρ is the number density of the molecules or atoms that form the liquid. B(r) is the bridge function. The analytic expression of B(r) is, in general, unknown. Unfortunately, the bridge function has to be approximated in practice. Though a number of theoretical and simulation procedures have been derived in recent years in order to get an approximate estimate of B(r) for different fluid models, the exact bridge function is not known for any system.  (Homeier et al., 1995)]. The solution of the OZ equations can be treated either as a fixed point problem for the unknown function γ(r), so that T OZ f (γ) = γ or as a multidimensional root finding problem F OZ r (γ) = T OZ f (γ) − γ = 0. In SASfit the fixed point problem T OZ f (γ) = γ is implemented as: As we assume that the system is isotropic, i.e. c(r) and γ(r) are only functions of the modulus of |r| = r, the Fourier transforms can be written as To calculate these transforms numerically on a discrete grid of points, the integration has to be substituted by a summation. Hereby, ∞ 0 dr is replaced by Np−1 j=0 ∆r, r j = (j + 1)∆r, q j = (j + 1)∆q, and ∞ 0 dq is replaced by Np−1 j=0 ∆q. The step width ∆r in real space has to be defined by the user as well as the number of grid points N p . The step width in reciprocal space is than obtained by ∆q = π (Np+1)∆r . It is important to mention that the values of the first elements of the arrays r i and k i equal ∆r and ∆q, respectively (r 0 = ∆r, q 0 = ∆q): c i = 4π∆r 2 (i + 1)∆q Np−1 j=0 c j (j + 1) sin π(j + 1)(i + 1) N p + 1 (77) γ i = 1 (2π) 3 4π∆q 2 (i + 1)∆r Np−1 j=0γ j (j + 1) sin π(j + 1)(i + 1) N p + 1 For calculating the Fourier transformation the library FFTW for fast Fourier transformation (Frigo & Johnson, 2005;Frigo & Johnson, 1997) has been used (source code available at http://www.fftw.org). The library supplies a discrete sin-transformation (type-I DST) named FFTW_RODFT00 which is doing the transformation: Np−1 j=0 X j sin π(j + 1)(k + 1) N p + 1 The unnormalized inverse of FFTW_RODFT00 is FFTW_RODFT00 itself. The same routine can be used for both transformations.
For the fixed point problem T OZ f (γ) = γ several iteration schemes are implemented.
Tests have shown, that no prediction can be made which one converges fastest or in case of e.g. very high volume or strong attraction, converges at all. Next to some simple iteration schemes also a the very efficient Anderson Acceleration for fixedpoint iteration (Anderson, 1965;Walker & Ni, 2011;Toth & Kelley, 2015) has been implemented. Besides these algorithms also several very efficient nonlinear system solvers based on Newton-Krylov solver technology, like GMRES (Generalized Minimal RESidual) (Saad & Schultz, 1986;Kelley, 2003), FGMRES (Flexible Generalized Minimal RESidual) (Saad, 1993), Bi-CGStab (Bi-Conjugate Gradient Stabilized) (van der Vorst, 1992;Kelley, 2003), and TFQMR (Transpose-Free Quasi-Minimal Residual) (Freund, 1993;Kelley, 2003) are supplied. These nonlinear system solvers are made available through the 'kinsol' library which is part of the SUNDIALS package (Hindmarsh et al., 2005) available at https://computing.llnl.gov/projects/ sundials. Test have shown that either the Anderson acceleration, GMRES or FGM-RES seems to be the most efficient algorithm to solve the Ornstein-Zernike equations with a minimum number of calls to the procedure evaluating the self consistent operator F OZ r (γ) or T OZ f (γ). The best choice for the iteration method depends often on the problem, i.e. the type of closure and the potential and also on the volume fraction.
Therefore, all methods mentioned above are still available in the program which might change after getting some more experience with them.

Parameters of numerical algorithms
Many numerical algorithms and their parameters in SASfit can be adjusted in detail if the particular application requires it. For example, the more efficient and flexible routine sasfit cubature for numerical integrations in multiple dimensions and a routine for optimized spherical averages (sasfit orient avg) was implemented. For both routines, an algorithm can be chosen by the user in a dedicated menu shown in figure 16 along with the respective parameters. Each integration strategy is also supplied as a multiple integration routine. Several of them are just nested calls and only few are designed for multidimensional integration. Because the menu interface should reflect this better, its revision is planned for the future. Fig. 16. A dedicated menu for choosing a numerical integration routine with its respective parameters especially for multi-dimensional integrations (sasfit cubature routine) and spherical averages (sasfit orient avg routine). This menu can be found in the fit or simulation window menu bar under 'Options' → 'Customize'.